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ABSTRACT 


This research considers a cantilever beam which can move 


axially in and out of a rigid frictionless hole and is free to 


vibrate laterally outside the hole. Two Euler equations 


describing the lateral and axial motion of the beam are 


presented. A transformation of coordinates to eliminate the 


moving boundary, and spatial non dimensionalization are used 
to transform the problem into a system of two coupled non 
linear partial differential equations with a fixed domain. A 


finite element formulation provides a numerical solution to 


the problem. Results for some problems are presented. 
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I. INTRODUCTION 


A literature search of the engineering journals shows that 
an investigation of the transient behavior of a cantilever 
beam, free to move axially in a frictionless hole at its 
‘fixed’ end, has not been undertaken to this date. In 1979, 
Boresi and Salinas prepared a report for the Naval Sea Systems 
Command, that formulates the problem and proposes a solution 
procedure. The report was the result of an interest in the 
transient behavior of a gun barrel during recoil following 
firing. [Ref.1] 

Hamilton’s principle was used to generate the governing 
partial differential equations for axial and lateral motion of 
the beam [(Ref.1]. As a result of axial motion of the beam, 
the length of the beam changes with time. Thus the ‘free’ end 
of the cantilever beam is a moving boundary point. If the 
beam is subjected to an axial force, then the beam length, 
that is the location of the ‘free’ end, is an unknown which is 
itself a solution of the problem. This is a ‘conjugate’ 
problem, wherein the boundary condition is a solution of a 
problem which can not be solved until the boundary extent is 
known. The analogy is of a dog chasing its own tail, or the 
‘catch 22’ syndrome. The dilemma is resolved by introducing 
a coordinate transformation which produces a classical two- 


point (fixed) boundary domain. The removal of the moving 











boundary is not without expense, as the resulting governing 


partial differential equations are significantly more 
complicated. Thus the complication of the boundary condition 
has been ‘transferred’ into the interior domain of the 
problem. The two equations governing axial and lateral 
motion, for beam length and lateral motion, are both coupled 
and nonlinear if the axial motion is not prescribed. 

Using the finite element method over the spatial domain, 
the two partial differential equations in space and time, are 
reduced to a system of ordinary differential equations in time 
only. That is, the original initial-(two-point) boundary 
value problem is transformed into a system of initial-value 
problems for the transient behavior at discrete points of the 
system. These nonlinear 0.D.E.’s are linearized using the 
quasi-linearization technique of Bellman [Ref.2], and then 
solved by using a fifth order Gear’ method for stiff equations 

This investigation adds further to the formulation of the 
problem by tne introduction of non dimensional variables. 
Additionally, the work also provides mathematical development 
and details required for the numerical solution of the 
problem. Restrictions and a generalization of the problem are 
also discussed. 

The scope of the problem suggests a cautious two-stage 
investigation. In the first stage, the axial motion as a 


function of time is prescribed. The result is the elimination 


of the need to solve the equation for axial motion. However, 








the equation can be used to solve for the axial force 
directly. Moreover, the remaining governing equation for 
lateral transient behavior is linear since the ‘length’ term 
in the equation is known. It is felt that the first stage 
investigation, which is the body of this thesis activity, 
would provide useful insight into the nature of the problem 
prior to undertaking the second stage investigation. In the 
second stage investigation, instead of prescribing the axial 
motion, the axial force at the sliding end is prescribed. As 
a result, the equation for the transient axial response needs 
to be solved in conjunction with the equation for transient 
lateral response, since now the length of the beam is also 
unknown. The second stage problem is formulated but not 


solved here. 


II. PROBLEM FORMULATION 


Consider the transient behavior of a cantilever beam 
fitted snugly into a frictionless hole as shown in Figure 
(2.1). The beam is free to move axially and laterally when an 
axial force F(t) is appiied, or when otherwise an axial 
displacement is imposed. The beam’s motion can be defined 
completely by its axial motion u(t) as a function of time, and 
its lateral motion v(x,t) as a function of both time and 
position along the x axis. Because of inertia, under certain 
conditions, such as when the axial force F(t) is a large 
magnitude impulse, the axial movement of the beam may tend to 
bend the beam by beam-column action or compress the beam 
axially by beam-bar action. These axial deformation effects 
are not considered here, that is u’ = 0. Therefore, it is 
assumed that all points along the x-axis of the beam 
experience the same axial motion. Thus, the instantaneous 
length of the beam, L(t), serves to describe the axial motion 
of all points of the beam. 

As the beam moves axially, the length of the beam outside 
the hole at any time t is defined as L(t). Because L(t) is 
changing with time the extent of the domain of the problem 
changes with time. It is this changing domain that results in 


the coupling of the equations which describe the lateral and 





axial motion of the beam. The changing domain is the essence 
of the problem and will be discussed at length in the 
development that follows. This investigation will be 
restricted to long slender beams, which in this case will be 
beams for which the length is equal to or greater than ten 
times either of the cross-sectional dimensions. With this 
restriction imposed, the Timoshenko Beam shear effects and 
rotary inertia, are neglected [Ref. 3}. However, as the beam 
length becomes shorter these effects become larger and loss of 


accuracy in the solution is expected. 





Figure 2.1 Cantilever Beam Moving 
Axially in a Frictionless Hole 














A. THE EULER EQUATION OF LATERAL MOTION 
Imposing equilibrium in the lateral direction and using 
small displacement theory results in the Euler Equation for 


the lateral motion of a beam, 


EI Vy45,(%,€)+ P V..(x,€) = p(x, t) 
t>0 (1) 


0< x < L(t) 


where the subscripts t and x denote partial differentiation 


with respect to time and position, respectively and; 


v(x,t) = the lateral displacement as a function of x and 
t. 
E = Young’s modulus of elasticity of the beam. 
I = moment of inertia of the beam cross-section. 
= the mass of the beam per unit length (constant). 
P = the internally applied load per unit length. 


The fourth order Euler Equation has two essential (forced) 
boundary conditions on displacement and slope at the ‘fixed’ 


left end, 


v(0,t)= 0 
(2) 
v,(0,t)= 0 











and two natural boundary conditions on moment and shear force 


at the ‘free’ right end, 


EI v,,(L,t) =M 
(3) 


where M and P are the applied moment and load, respectively. 
The homogeneous boundary conditions (M = P = 0) are the 
boundary conditions considered here. However, a verification 
of the solution method is presented where the non homogeneous 
boundary conditions are imposed. The term ‘fixed’ end is used 
in reference to the boundary located at the left end of the 
beam’s domain (See Figure 2.2), i.e., at x=0. As a result of 
the axial motion, the point on the beam at this left or 
‘fixed’ end is changing with time. 

The natural boundary conditions at the free end (Eqs. 3), 
for moment and shear, occurs at the right end point of the 
beam (i.e., at x = L) for all time t. It is the fact that the 
argument L in Equations (3) is changing with time that makes 
the natural boundary conditions troublesome. These so called 
moving boundary conditions (or changing domain) will be 
discussed later at length. 

The Euler Equation for lateral motion is also a second 
order differential equation in time. To obtain a solution, 
two initial conditions, one on its lateral position v(x,0), 


and one on its velocity U,(x,0), along the x axis will be 








needed. These initial conditions will depend on the specific 


problem being solved. 


B. THE EULER EQUATION OF AXIAL MOTION 

If F(t) rather than L(t) is prescribed, then a 
differential equation defining L(t) is needed. Again, using 
principles of equilibrium for motion in the axial direction, 


the following Euler equation for axial behavior is obtained, 


1 


F(t) 4 
or (4) 





3s 1 2 = 2 = 


Equation (4) is subject to the initial conditions, 


L(0) = L, 
(5) 
L(0)= £, 


where L, is the initial length of the beam at time t = 0. The 
dot and double dot above L denote the first and second 
derivatives with respect to time respectively, that is the 
velocity and acceleration of axial motion. 

Together Equations (1) and (4) along with their respective 
boundary and initial conditions form a coupled and nonlinear, 
initial-boundary value problem. When the force F(t) is known, 
these coupled nonlinear equations can be solved using the 


finite element method with a linearization scheme to find 











u(x,t) and L(t). When L(t) is specified, Equation (4) yields 
F(t) directly. 


C. THE MOVING BOUNDARY 

In the boundary conditions described in Equations (3), the 
beam length L(t) is a function of time. Thus the boundary 
conditions are conditions on a boundary which is moving. 
Graphically this is shown in Figure (2.2). The curved 
boundary of the region of integration presents a problem. The 
desire is to remove the argument of time varying length from 
the boundary condition at the free end. In essence, we desire 
to secure the boundary. Graphically the boundary becomes a 
Straight line where previously it was a curved line (See 
Figure 2.3). This can be achieved by using a coordinate 
transformation as shown in the next section. 


t 
x=L(t) 






Vy = Vxx = 0 


V= F(x) 


Figure 2.2 Region of Integration for Equation (1) 












Ves *Ve55=0 


v= F, (5) I. 


Figure 2.3 Region of Integration for Equation (25) 


D. THE COORDINATE TRANSFORMATION 


New variables (x,t) and N(t) are introduced as follows: 





a L(t) (6) 


It should be noted that § is a non-dimensional variable with 


respect to the spatial domain. The lateral deflection now 


becomes a function of these variables as shown below. 


v{ x(6,n),t(6,n)) , Or (7) 
v{ & (x,t) -N(x, t)) 








Considering the relations defined in Equations (6), the 


following partial derivatives are obtained. 


ala 
Ox L 
5 rxL - $b where L - OL 
ot Ee ot 
8 
an oo (8) 
ot 
M = 9 
Ox 


1. Transformation of the Spatial Fourth Derivative of v 
Considering Equation (7), the transformation of the 
spatial fourth derivative on lateral displacement, v,,,, to the 
new coordinate § is accomplished through a _ series of 
differentiations using the chain rule. The first 


differentiation results in, 


dv _ av 0& 


Ox of dx om ox (9) 


dx 0&6 dx an ox 


Following the substitution of Equations (8) into Equation (9), 


Equation (10) is obtained. 


(10) 


11 





After another differentiation with the chain rule the second 


spacial derivative is found. 





- selz" ee . lz Ee (11) 


Again, using Equations (8), the second derivative is equal to, 


vy =tyv (12) 


Likewise, the third derivative is, 


a A 
Vax ~ Ta Vane (23) 


and finally the fourth derivative is, 


SA 
xox a Verte tt4) 


2. Transformation of the Time Second Derivative of v 
The transformation of the time second derivative on 
lateral deflection (or acceleration), UW. to the new 
coordinates § and n is performed in a similar fashion as was 


the transformation of it spacial derivatives. Once again, 











using Equation (7) and the chain rule, the following 


expression for the first time derivative is obtained. 


06 ot oandat (15) 


Substituting the Equation (8) values of the partial 
derivatives into Equation (15) results in the following 


expression for V, 


Vv. = _ SL, + Vv (16) 


Another time derivative using the chain rule results in the 


following equation for Vy, 


(v. } 


ein 
| Sy, “i (17) 


ety, + Vi, a 2. Ey sy, [on 


tt 


Es 
ot 
ao 
ot 
= 
a& 


Again, using Equations (8), 


Obie EE yea. Sons 
de EY Oe ee 


13 








along with the product rule of differentiation, we obtain 


(19) 


Recalling that the coordinate transformation on time stated 


that t=, it follows that, 


. : OL oL 
L(t) = Z£(n) = LE = 22 
De ae On (20) 


Now, using the quotient rule of differentiation, Equation (19) 


becomes, 


(21) 


Finally, after multiplying and collecting like terms, v, 


becomes, 


14 








E. THE FINAL EULER EQUATIONS 
Using the transformed operators, the Euler Equations are 
rewritten in terms of the new coordinates, & and 1. 
1. The Transformed Euler Equation for Lateral Deflection 
Substituting the transformed operators from Equations 
(14) and (22), into the original Euler Equation for lateral 


deflection, 


ElV pcx + PYee = P(X, €) (1) 


results in the following Euler Equation transformed to the & 


and 1) coordinates, 


EI 2 |Z L ig tle _¢i i, 
Fa Very * P E EE Jv + 2§ z}y 2b 5 Ven 65%; + Van |= 


15 








Multiplying through by the inverse of the coefficient of Ve: 


gives, 


pL* L + L - L - L + = 
Yeu * er [* E Jv 26 E Mee Na ve 
xf 25 


0<&<1 


o<yn 


and its boundary conditions, 


v(0,n) = 0 Vig (1,1) = 0 


(26) 
v, (0,1) 5 0 Vere (1,1) oe 0 


The boundary conditions are now functions of & over the domain 
0< €& <1, in lieu of x over the domain 0 < x < L(t). The 
initial conditions on deflection and velocity will be 
functions of € as well. 
2. The Transformed Euler Equation for Axial Motion 
The coordinate transformation on the Euler Equation of 


axial motion shown again here, 


1 


tt) * Sor, 





[EIvi. (L,t) - pve(L,t)] = — F(t) (4) 





16 





results in the transformed equation, 





2 
», 1 |[BIy2 en ae = 
L pL, ser Veg tee of Evan + yqc1m | 565 
1 
F(n) 
pr, 


subjected to the initial conditions in Equations (5). 
3. Non dimensionalization of the Lateral Deflection, v 

The purpose of the coordinate transformation just 
completed was to deal with the difficulty presented by the 
moving boundary condition at the free end of the beam. The 
four boundary conditions of Equations (2) and (3) were also 
transformed to the € and nN coordinates as shown in Equations 
(26). One of the great difficulties encountered in this 
investigation resulted from the coordinate transformation 
performed on the boundary conditions. After the introduction 
of the non dimensional variable &, the finite element method 
(FEM) of Chapter III was pursued. This included an attempt to 
confirm the FEM program on a couple of statics problems with 
known solutions. The resulting FEM solutions were L? larger 
for displacements and L? larger for slopes. We had simply 
imposed the load (or moment) as one would have if the problem 
had the dimensional independent variable x, when in fact x had 


been replaced by the non dimensional variable § = x/L. This 


17 





problem was eventually resolved by introducing the non 


dimensional displacement, v*, defined as, 
(29) 


Its derivatives, 


ble 


(30) 


are used such that, 


= 9 « 9 (rye) = LV; 


Vv as 
5 ECE, 
, ; (31) 
Vee = aa = eo )= LViE 


and in the same fashion, 
(32) 


= Lie 


and, 
(33) 


Veeee = LViee 
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After making the substitutions into the equations of 
‘lateral and axial motion, Equations (25) and (28) 
respectively, the spatially non dimensional Euler equations 
are obtained. 

@. Final Euler Equation of Lateral Motion 

. Es Ps . os . 
Vizgs * B E Ss Vent 26 Ve - 26LViq - EL VE + L Vig |- 


0<E<1 
0o<7 


where, 


(35) 
eee ng 
P =I p(§,n) 
and its boundary conditions are, 
v°(0,n) = 0 Veg(1,n) = 0 
(36) 
vz (0,1) = 0 Vere (1,N) = 0 


Again, discussion of the initial conditions will be 


delayed until later. 
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b. Final Euler Equation of Axial Motion 





1 (3 


= vee(1,9) ~ p[-Lvg(1,q) + ZL va(1,n)]’ z 


Z. 
pr, 7 (37) 


o<&<1 
0o<7n 


and initial conditions, 


L(0) 


tt 
i 
° 


(38) 
£(0) 


" 
tH. 
° 


F. CASES 

There are two general cases for which the transient 
behavior of the cantilever beam may be considered. Recall 
that the beam is free to move axially when an axial force F(t) 
is applied resulting in an axial displacement, or when an 
axial displacement L(t) is otherwise imposed. It was shown in 
the Euler equation of axial motion (Eq. 37) that the axial 
motion described by L(t) depends on F(t). However, if this 
axial motion is specified simply as some function of time 


alone then the problem is greatly simplified. 
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1. Case One, L(t) Prescribed 

If L(t) is known then the problem is reduced to 
finding a solution to the Euler equation of lateral motion 
(Eq. 34) subject to its boundary and initial conditions. The 
problem is a linear, initial-boundary value problem. 

An even further simplification occurs if the axial 
motion is so slow that the time derivatives of L(t) are 
negligible. Certain linear functions of L(t) can conveniently 
provide such a condition where the rate (or velocity) is made 
small, and the second time derivative (acceleration) is always 


zero. In this case Equation (34) reduces to, 


+ 


* pLi oo. _ 
Vecee Er Ym = 0 139) 





subject to its boundary and initial conditions and where 
p*(&,n) = 0 for a beam with no internal loading. Note that L‘ 
is not a constant here, since L(t) is a prescribed function of 
time which needs to be known. 

In either of the cases (Eqs. (35) or (39)), a closed 
form solution to the equation(s) is not possible. The finite 
element method (FEM), to be presented in Chapter III, was used 
to obtain approximate solutions for both of these cases. 

2. Case Two, F(t) Prescribed 
In the case where F(t) is prescribed, both Equations 


(34) and (37), subject to their respective boundary and 
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initial conditions must be solved. We recall, that together 
these two equations form a coupled, initial-boundary value 
problem. Moreover, they are both now nonlinear, as they 
contain terms with both dependent variables, L and Uv and their 
derivatives. Therefore, it is necessary to linearize both 
equations. 

Any number of different strategies are possible for 
the linearization of these equations. The strategy used here 
will be to treat each dependent variable as a ‘primary’ or 
‘secondary’ variable in accordance with the following scheme. 
The assignment of ‘primary’ or ‘secondary’ status will differ 
depending upon which equation is to be linearized. In the 
linearization of the equation of lateral motion, the ‘primary’ 
variables are lateral deflection, v and its derivatives; and 
the ‘secondary’ variables are axial motion, L and its 
derivatives. As ‘secondary’ variables in the equation for 
lateral motion, L and its derivatives are evaluated at the 
previous time. On the other hand, for the equation of axial 
motion, L and its derivatives are considered the ‘primary’ 
variables and UV and its derivatives are the ‘secondary’ 
variables. In this case VU and its derivatives are evaluated 
at the previous time step. 

The linearization of Equation (34) is quite simple 
because in the nonlinear product terms, the primary variables 
(v and its derivatives) appear linearly. Therefore, it only 


becomes necessary to evaluate the secondary variables in these 


22 





product terms at the rrevious time. For completeness the 


linearized equation is shown below, 


2 


2 
Views + BL 1S? | Vie + 26(2 | vi - 262, vin - GL, Ve + L, Vag |= 


p’(&,n). 


o<E<1 


0o<7n 
(40) 


where the * subscript on a variable (or term) denotes that 
the variable (or term) is evaluated at the previous time and 
therefore is not an unknown in the equation. 

The linearization of Equation (37) is not so simple 
because the primary variables (L and its derivatives) do not 
appear linearly in the nonlinear product terms. 


If Equation (37) is expanded, 


2 
~, 1 [vee | oa pee yee 1 - ee 
7 Dp, ee | zz, | ve*(1/)] + Zz [LL ve* (1/7) * vad -N) | 


1 2 v2 = 1 
- dy, [z? ve7(1,0)] = F(n) 
3 [ y (1,7 | Zpl, | (aa 


each of the bracketed [ ] terms are nonlinear. These terms 
can be linearized in a number of different ways. Since this 


is primarily an "L" equation, the "L" operators will be 
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linearized using the quasilinearization technique or Bellman 


& Kalaba [Ref.2]. The nonlinear terms then become, 
Vee (1,1) a 2 
1 = Vise (= ~ 3") 


where, Vic(1,n) = |--© v.(1,n) + 4 @,(1,n) : 
’ i ’ i . f 2m . t (42) 
2. 2? vi7(1,n) =[0.(1,9)]? (-£2 + 22, 2) 
3. LE vi(1,N) va(1,-q) = 8(1,N) V.(1,N) LL 


4. L? vi(1,n) = v2(1,n) (-L? + 2b, L) 


where 1, is the length of a finite element after the beam is 
discretized, 90 represents slope v;,. Recalling that n=t, in 
these equations, the subscript nN denotes partial 
differentiation with respect to time. Again, the * subscript 


on a variable (or term) denotes that the value of the variable 


(or term) from the previous integration is used. 





Equation (37) can be rewritten in terms of the L operators 


‘and the following groups of terms, 








_ EI. 3 [_ 6 4 f 
C, = pL, x Zz v.(1,7) + em cm | 
2: HE sce 2 4 ‘ 
C, a ZpL, ra Tr v.(1,7N) + I. ecm) | 
2 1a 2 2 
C; 2L, 6.(1,7N) Li 
C, = -—4 02(1,n) 2, 
2L, (43) 
C= = O.(1,0) V1) 2, 
0 
sy 1) ae og igyee 2 
C a ve(1,N) L 
Pony ee 
C, os Vv. (1,7) q, 
= 1 
ae 
Using Equations (43), Equation (37) becomes, 
Dt C+ CL: Cy +0, Lb * Ch. 4-C. 4 Che =O F (y) (44) 


Equation (40), its boundary and initial conditions, 
and Equation (44) with its initial conditions, now form a 
coupled, linear initial-boundary value problem. A closed form 
solution of these equations is not possible. The finite 
element method (FEM), to be presented in Chapter III, could be 


used to obtain approximate solutions to these equations. 
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III. SOLUTION METHOD 


A. FINITE ELEMENT METHOD DEVELOPMENT 

Considering the L(t) specified case first, the task is to 
solve Equation (34) together with its boundary conditions 
given by Equation (36), and its initial conditions as 
determined by the problem being investigated. Definition of 
the initial conditions will require further discussion which 
will be conducted later in this development. 

Equation (34) is a linear partial differential equation in 


one unknown, Vr (E,t) when L(t) is specified. Recalling that 


N=t, here t will replace 17. An approximate numerical 
solution of this equation together with its initial and 
boundary conditions can be obtained by a Galerkin finite 
element formulation. 
1. Construction of the Beam Elament 
The fourth order system of Equation (34) requires C’ 
continuity (continuity of function and slope). In order to 


obtain C’ continuity, an element with deflection v*, and slope 
6°, (8° will represent v," in the development that follows, as 


degrees of freedom (DOF) at each end point is required. This 


means each element must have a minimum of four DOF, which 


requires a cubic polynomial. These interpolation polynomials 








are the set of cubic spline shape functions listed below and 


detailed in Appendix A. 


Lp a, Sed ge a 
Be 

gp eee ge as 
ao ee ae 

(45) 

3 42 2 43 
q, = = eo 
ae a 

3). Be 1 
ee ae 


Figure (A.1) shows that shape functions q, and q,;, are 


associated with the displacement DOF (Vv; ,V; where subscripts 


1 and 2 represent node points 1 and 2) at the element end 


points; and the even numbered shape function q, and q, are 


associated with the slope DOF (8; , 62) at the same locations. 


2. Global FEM Formulation 


In terms of global shape functions Q,, the FEM 


approximation v* to 9“ is given by, 


N 
vi=et=Q975 =) Q, 5, (46) 
izl 


where N is the number of elements, N = 2N +2 is the number of 
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global DOF, and 6, are the global degrees of freedom. The 


global degrees of freedom, 5, are defined as follows, 


87 = (8, 8, : 8, 8, : 8 8, 2 ~ : 8s 85) (47) 


where subscripts 1,2,3,...(N-1) are DOF identifiers. In terms 


of displacements and slopes, we have, 


67= = (v, 6, i V2 8, 7 V3 8, bo 3 Vue Oy.2 ) (78) 


where the subscripts 1,2,3,...(N+t1) refer to the Global Nodal 
Points (GNP) and, < >* indicates the non dimensional forms of 
v and 6, that is v* and @*. 


The relationship between the global degrees of freedom 
5, and, 0° and 6° ; is such that for odd i (1,3,5,...N-1); 
6, = 0° at GNP [(it+1)/2]. (i.e., 8, = V°(GNP 1), 8, = 9° (GNP 
2), «.-) For even i (2,4,6,...N); 5, = 6° at GNP [i/2]. 


(i.e., 5, = 6° (Gnp 1), 8, = 6° (GNP 2), ...) 


Each i** GNP has two global shape functions, and hence 


two DOF. An odd numbered shape function Q,;, gives 


displacement V* at GNP{ (i+1)/2} = Vicierys2) = 5,, and an even 
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numbered shape function Oi; gives slope 6: at 


GNP{i/2} = 9j4;2) = 5,- 


3. Galerkin FEM Formulation 
In accordance with the Galerkin Finite Element Method 


(FEM), we form the approximate solution of vU, 


v°(&,t) = 0°(E,t) = O7(&) 5(t) (49) 


Using the above approximation, the residual function for the 


Euler equation of lateral motion (Eq. 34) is, 


R(§,t) = 2(9"} - p* (50) 


or, 


R(E,t) = (9)... + sons 


oe ac ae ba re Pee. 2 Be ee ante |oes 
| = Yee 26-9 2EL (Vv F EL (Y ), + LV | P 
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After the final substitution, the residual is, 


R(€,t) = (97 _ 5 + 


Pat T 2 L? T T r” T T 
pH lor), 8 + 26% (or) 8 - 262 (or), 5 - gz lo"), 8 + 20° | 


- p* 
(52) 


The Galerkin finite element equation is obtained by requiring 
that the residual function be orthogonal to each of the basis 


functions. That is, 


ff or ab = 0 (53) 
0 
Substitution of Equation (51) into Equation (52) gives, 


e T + a - T + 
f, ele leeee A 3 + BS j, & oo"), ab 8 
(54) 


2p {* bo (or), a& 5 - 22 [* Bolo"), ak 8 - 
Bi {° Eo(o"), ab 8+ Bx f olor) a& 5 = [op ak 
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After performing two integrations by parts on the first term 


(See Appendix D), Equation (53) becomes, 


B.T. +f 2, (o7),, ab b+ pe f, E2 2 (07),, dé 5 + 
je f 60 (07), a& 5 - 2B f- elo"), db (55) 


Bi f° Bolo"), d& 5+ Bx [' o(o") dS = [opr ak 


where B.T. is the vector of boundary terms resulting from the 


two integrations by parts, 


B.T. = Q(%") Lo - 219"), [0 (56) 


The kronecker delta property of the shape functions Q,, 


results in the reduction of the boundary term vector to, 


— (97) 55, (0) 
(9"),, (0) 
0 
0 (57) 
0 
0 
(97) e— (2) 
=t9"),, (2) 
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The non zero terms in the B.T. vector represent the non 
dimensionalized loads and moments applied at the boundaries. 
If there is no applied load (or moment) at the boundary, then 
B.T. is simply the zero vector. 

A discussion of the second term of Equation (54) 
follows. The integration by parts performed on the first term 
of Equation (54) results in a self adjoint (symmetric) 
operator, a condition which is generally desired since it 
reduces storage requirements during computer processing. 


Integration by parts on the second term, 


[8 200% 4 a6 8 (58) 


results in a B.T. on 6° evaluated at &=1. Since the value 


of 6° (1) is unknown, an integration by parts was not performed 


on the Vee. 
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Letting, 
a= {'9,, (97), a& 
B= {' & olor), at 
c= [* Eo(o), a€ 
p= f alo"), at 


1 
r- {op d& - B.T. 


the final Galerkin Equation is, 


ad+B+ 55+ 2p cb - 


2prc5-Bic&+BPLDb =F 





The details of the construction and form of the A, B, 
C, and D matrices are contained in Appendix B. 
4. Integration of the Galerkin Equation 
a. Reduction of the Second Order System 
Equation (59) is a system of second order ordinary 
differential equation in time. For numerical integration 
purposes, it is desirable to reduce this second order system 


to a first order system in time. 
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For compactness, let, 





L 
uz 22 ¢ 
L 
41 
P al 


Then, Equation (59) becomes, 


Dd5=M5+K5+P 


Letting @ = 3, it follows that, 


and Equation (59) now becomes the first 


equations, 


D® = Md + K5+P 
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(61) 


(62) 


(63) 


order system of 


(64) 





In explicit matrix form, Equations (62) and (63) may be 


written, 
4 BY 
zr:0 $ Oo: 2 {id 0 (65) 
= + 
0 D hh) K M |i P 
Letting, 
r:0 0 =I 
Gt 4 = , A= i= 
0:oD Ki M 





(66) 


the second order system of Equation (59) is reduced to the 


following first order system in time. 


Gi=HA+Q wh 


Vector Q, becomes the zero vector if vector P is a 
zero vector. Referring to Equations (59) and (60) we see that 
P is actually defined by vector F which is defined further by 
the boundary term vector B.T., and the vector describing the 
contribution of an internally applied load, p*. If B.T. is 


the zero vector (no applied moment or load at the boundaries) 
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and there is no internally applied load, then P and Q are zero 
vectors. Evaluation of the B.T. vector has also resulted iu 
satisfying the natural boundary conditions imposed on Equation 


(34). Equation (66) now becomes, 


Gh=8X Ree) 


b. Boundary and Initial Conditions 

Prior to integrating the system described in 
Equation (67), the boundary and initial conditions on Equation 
(34) must be imposed. ‘The boundary conditions at the free end 
(E=1), were imposed through the boundary term vector as 
previously described. The strategy used to impose the 
essential boundary conditions at the "fixed end" (€=0), is one 
in which the deflection and slope at the "fixed end" are set 
to zero when the A vector is initialized, and the G and BH 
matrices are altered such that the deflection and slope at § 
= 0, remain constant with time. If the first and second time 
derivatives of deflection and slope at € = 0 are constant and 
equal to zero, then the desired conditions of zero deflection 
and slope at § = 0 are obtained providing that the initial 
conditions on deflection (0(0,0) = 0), and slope (8(0,0) = 0), 
are satisfied. 

Initial conditions are imposed through the 
initialization of the A vector in accordance with the problem 


being investigated. Referring back to the global FEM 
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formulation, we recall that each nodal point has two DOF. To 
‘satisfy the two DOF, deflection (and its velocity), as well as 
Slope (and its velocity), must be initially defined at each 


node. The initial conditicus mest alsu satisfy the essential 


boundary conditions at § = 0. The initial conditions are more 


clearly understuod if the A vector is given in greater detail, 


Vay (69) 


where the * subscript indicates the terms are the non 


dimensional variables (and their derivatives) and therefore, 


are functions of §, not x. 
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The matrices and vectors of Equation (67), modified 


for the boundary and initial conditions follow, 


(70) 


where only the first two rows of D are altered as shown and, 





(71) 


where only the first two rows of I, K, and M have been altered 


as shown and, 


aoe uo 
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and, 


(73) 
r 5 eee - 


B. FEM VERIFICATION 


The static cantilever beam provides a problem for which a 


known solution is available for comparison and verification of 


the FEM development and Fortran code. 


For the static cantilever beam, the Euler equation of 


lateral motion (Eq. 1) is reduced to, 


EIV,...(*) = p(x) O<x<L (74) 
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with the boundary conditions, 


v(0) =0 EIv,,(L) = M (moment) 
(75) 
v,(0) = 0 EIV,,,(L) = P (load) 


Referring to the Euler equation of lateral motion (Fq. 


25), and considering that for the static cantilever beam, 


Le =0 (76) 


Equation (73) transformed to the non dimensional coordinate §& 


becomes, 


4 
Vests = zp P (6) o<&<1 (77) 


with the boundary conditions, 


EI 


v(0) = 0 Fay (l) = M 
(78) 
v,(0) = 0 A Vag(l) =P 


Referring to Equation (34), the final spatially non 


dimensional static beam equation, where the lateral deflection 
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has also been non dimensionalized, becomes, 


3 
Vexee(§) = arP tb) O0<E<1 (79) 


with the boundary conditions, 


v"(0) = 0 =e vig (1) =M 
(80) 
vi(0) = 0 Fa Ving (1) =P 


The Galerkin FEM formulation for Equation (78) and its 
boundary conditions is obtained from Equation (59). Again, 


recalling that, 


L=L=0 (81) 
5 =0 
Equation (59) becomes, 
ASi=F (82) 
where, 
F = Qp'd& -B.T. (83) 
0 
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and, if there is no excitation internal to the system other 


than at the boundaries (p* = 0) then, 


F=0 -B.T. (84) 


The boundary conditions given in Equations (79), must be 
imposed prior to solving the system of Equation (81). To do 


this the boundary term vector B.T., resulting from the 


integration by parts on the Vie: operator and shown in 


Equation (56), 


0 (56) 
0 


Vege (1) 
-Vie (1) 


is evaluated using the boundary conditions at the free end of 


the beam (€ = 1). 
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Rearrangement of terms in Equations (79) gives, 


(86) 


The B.T. vector can be rewritten as, 


~Viee (0) 
Veg (0) 
0 


B.T. = ; (87) 


Next the -oundary conditions at the fixed end (& = 0) are 
imposed. Recalling that 6, = v; = v"(0) and 6, = 6; = @°(0), the 


boundary conditions at the fixed end are imposed by altering 
the first two rows of both the A matrix, and the B.T. vector 


to force 6, and 8, to zero. 
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The final system can be written as, 


100 000 5, 0 
010 00 0 5, 0 
Ay; = Aun 
= (88) 
PL? 
Ora TEE 
ML 
Ay, Aw 55 Er 
where N is the number of degrees of freedom. 
The solution to this system is, 
Vi 
ee al (89) 
Vy 
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PL? ML 


If “= (or 22) is set at unity, then to obtain the 6 for 
EI EI 
2 2 
actual values of = (or =) , the 6 vector is multiplied by == 
ML 
(or =a): 


The exact solutions for the deflection and slope at the 
free end of a cantilever beam subject to a concentrated load 


(or moment) are obtained from the following expressions, 


3 z 
(90) 
2 


where P (or M) is the load (or moment) applied at the free 
end. 

The 6 vector is the vector of non dimensional deflections, 
and slopes. The dimens 1al vector is obtain by multiplying 
the non dimensional displacements (6, for odd i) by L in 
accordance with v = Lv’. Since slope is a non dimensional 
quantity to begin with, the 6°’s (6, for even i) are equal to 
the 0's. 

The Fortran code used for this verification and comparison 
is located in Appendix C. Results of the verification, also 


in Appendix C, confirm the FEM development. 
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C. THE F(t) PRESCRIBED SOLUTION METHOD 

The final and most complex case posed in Chapter II was 
the case where the axial force, F(T) is prescribed. In this 
case an equation for axial motion, in addition to the Euler 
equation of lateral motion, was required. That equation was 
the Euler equation of axial motion. Together these two 
equations form a nonlinear, coupled, initial-boundary value 
problem. After the linearization of these equations, the 
linear, coupled, initial-boundary value problem consisted of 


the following equations, 


2 2 
Veet + B. 1G? | Viz + 26(2| VE - 2E2, Vin - Ef, Vi + L, Vay 


= p’(&,n) = 0 
(91) 
L+ce,+0,b+C,+¢C£+ C+ C+ CL = CF) (92) 
where p* = 0 for the no internal excitation case, and C, are 


defined in Equations (43). 
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By defining new terms, J, which require updating during 


the time integration process, 


J,(t) = C, + GO, + C, 
J,(t) =c, (93) 


J,(T) =C, + Cy + 


Equation (45) becomes, 


L+J,L+J,L2=C, F(t) - J; (94) 


Equation (93) is a second order differential equation in 


time. Letting, 


(95) 


and using substitution, the following system of two first 


order differential equations is obtain; 


Z=L (96) 


P+, L+I,L=C, F(t) - J, 


Since L(t) and its derivatives do not depend on a spatial 
variable, Equations (95) do not require a FEM formulation and 


are directly added to the system of equations for lateral 
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motion. The matrix equations for lateral motion will be 
‘similar in structure to the system of Equation (64). However, 
the ‘secondary’ variables (or terms) introduced in the 
linearization of the equation, are assigned their values from 
the previous integration. Therefore, the sub matrices D, K, 
and M which contain these variables (or terms), appear with 
the * subscript. Q is the zero vector because there is no 
internal excitation, and the natural boundary conditions used 
to evaluate the B.T. vector are equal to zero (moment=load=0). 
The matrix system of equations which is obtained after adding 


Equations (95) to the system in Equation (96) is shown here, 





The final system for the F(t) prescribed case is, 


0: 0 
I 0 § 
- >. 7 
0 Dd, rt) = 
0: O 
0 Sed 0 1 0 L 
0 oS 0 :0O0: i121 & 
0 0 ‘ 
0 P| é 0 
if t 
(99) 
K, M, @ + 0 
0 0 
0 2 ¢) 0 1 L 0 
0 “ 0 -J, =J5 g¢ C,F - Jd, 


The boundary and initial conditions for the equation of 


lateral motion, and the initial conditions for the equation of 








axial motion, are imposed using the strategy applied in the 
L(t) specified case, and described in Section A of this 


Cha»ter. 
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Iv. CASE STUDY REPORT AND CONCLUSIONS 


A. GENERAL DISCUSSION 

After verification of the finite element code the 
investigation of the transient problem began. The primary 
emphasis of the studies that follow is on obtaining solutions 
to problems, and not on investigation of numerical 
considerations. However, when appropriate the researcher's 
thoughts on such considerations will be presented. 

The case studies reported are investigations of the L(t) 
prescribed condition. For reference, the system developed in 


Chapter III for the L(t) prescribed case is repeated here, 


r:o {is Oo: x |{d (100) 
0: DI Ki: Miw 
or, 
Gi = BA e103) 


The above system does not reflect the alterations made to 
impose the boundary conditions on load and moment at the free 


end as the case studies addressed only the case of homogeneous 
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boundary conditions (P=M=0), with no internal excitation 
(p’=0). Thus the P vector is the zero vector, and does not 
appear in the system above. 

The transient problem introduces the requirement for a 
numerical integration method. To perform the integration on 
the system above, the IMSL, Inc. integration subroutine, 
DIVPAG was chosen. DIVPAG is a double precision first-order, 
initial-value, ordinary differential equation solver. 

Two classes of implicit linear multi step methods are 
available. The first is the Adams’s method and the second is 
the backward differentiation formula (up to fifth order), also 
called Gear’s stiff method. An accepted measure of stiffness 
is the ratio of the maximum and minimum eigenvalues (Aya, Ania) 
of a system. A problem is considered stiff for very large 
Anox/ Amin Yatios. The vibrating cantilever beam equation of 
motion results in a stiff system, and therefore, Gear’s stiff 
method is used. 

1. Time Step Convergence 

The integration routine uses an internally determined 
time step such that a measure of global error does not exceed 
a user specified tolerance. This feature of DIVPAG provides 
error control to the user of the integration routine. 
However, a recognized short coming of this integration package 
as applied to this problem, is the inability to update the G 


matrix on the left hand side of Equation (100) at each of the 
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subroutine determined time steps. That is, as a result of 
the call structure of the subroutine, an update of the finite 
element matrix G for a change in L(t) and its derivatives, is 
only possible outside of the subroutine. For this reason 
DIVPAG is placed in a Do Loop and the G matrix is updated at 
each entry to DIVPAG. Although the H matrix could be 
evaluated inside DIVPAG via a FNC subroutine argument of 
DIVPAG, in this investigation it was not. It was updated at 
the same time the G matrix was, that is, at each entry to 
DIVPAG. The accuracy of the numerical solution depends upon 
the frequency of updating of the G and H matrices. Entries to 
DIVPAG were at .025 second intervals for all case studies with 
the exception of the final case study for which entries were 
made at .01 second intervals. A rapidly changing L(t) 
requires more updating of the matrices than would a slowly 
changing L(t). In effect, a solution ultimately should be 
checked for "time grid" independence. 
2. Spatial Grid Convergence 

Convergence of the solution for the spatial grid is 
yet another consideration. The solution is a function of the 
number of elements (i.e., NDOF). For linear problems, it can 


be shown that in the limit, as the number of DOF approaches 
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infinity, the approximate solution of V approaches the exact 


solution v, that is, 


LIM ¥ (102) 


However, for a nonlinear problem there is no guarantee, but a 
likelihood that the approximate solution converges to the 
exact solution in the limit as the number DOF approaches 
infinity. A preliminary study conducted during the first case 
study showed that negligible differences existed between the 
eight and sixteen element solutions for that particular 
problem. This was the basis for the use of an eight element 
solution for all subsequent problems. However, it is 
recognized that the FEM model changes with length (or time). 
Since the number of elements (NEL) is constant with time, 
convergence for a given NEL may change with time as well. 
Furthermore, the effects of material properties, geometric 
dimensions, and functions of L(t) (and its derivatives), may 
also influence the NEL required for a spatial, grid 
independent solution. 
3. Computational Effort 

Related to the stiff character of the problem, is the 
very large amount of computational effort (CPU time) required 
to obtaia solutions. A study of CPU requirements was not 


conducted. However, integration of a problem over a real time 
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five second period took as long as a week. Typically, DIVPAG 
performed its integration over 2.‘*® second steps. Thus, every 
-1 second increment in time required approximately 500,000 
integration steps. Processing was conducted using the Naval 
Postgraduate School IBM 3033 main frame time share system 
during weekday non peak hours (1800-0900), and weekends. 
During these periods, it is estimated that approximately 20 
percent of time share CPU was allocated to the processing of 
this job. A restart capability was coded to assist in 
processing during non peak hours only. 
4. The Case Study Beam 

The case studies that follow are conducted using 
material properties similar to those of plexiglass. The 
modulus of elasticity (E) is equal to 100,000 psi. The 
initial length of the beam, L, is 10.0 inches for all case 
studies. Two moments of inertia of the beam cross-section 
(indication of beam rigidity, which effects the stiffness of 
the problem) are obtained from the cross-sectional dimensions 
shown in Figure (4.1). In recognition of the stiff nature of 
the problem, and in interest of solving a realistic problem in 
the minimum amount of time possible, our desire was to select 
a material with the smallest value of frequency, which is 


proportional to, 


EI 


ryh (103) 
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Thus, for a beam of given geometry (I/L‘), plexiglass was 


selected as a realistic material with the smallest E/p ratio. 


0.125" 
0.0625" 





Figure 4.1 Beam Geometry 


B. FIRST CASE STUDY, NEGLIGIBLE DERIVATIVES OF L(t) 

The strategy used in the case studies is to progress from 
the less to the more ‘difficult’ cases. What is intended by 
‘difficult’, is that fewer differential, and hence, finite 
element operators are involved in the Fortran coding for the 
less difficult case. The general program development logic is 
the same for all the transient case studies, however; it is 
generally good practice to limit the size of the code until 
the logic is tested and functioning as expected. By 
eliminating the derivatives of L(t), only the A part of the K 


matrix (See Equation (105)) and the D matrix of Equation (99) 


need to be evaluated. 
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If the derivatives of the specified function of L(t) are 


zero or negligible, the equation of lateral motion is, 


: pLt .. _ 
Vous * gr Van = ° e104) 





subject to boundary and initial conditions. 


For this case, the matrix system of Equation (99) becomes, 


(105) 








ee om (22 ; tJe| (106) 


is reduced to, 


K=- A (107) 


This case was examined for the plexiglass beam with the 


larger cross sectional dimensions. The material and geometric 


properties used are p(mass/unit length) = 6.988E-6 lbf-S?*/in’ 


= 81.38E-6 in‘. 


(slugs/in) and, moment of inertia I 








The homogenous boundary conditions are imposed as 


described in Chapter III. The initial conditions (A vector) 
are imposed by an initial parabolic deformation (Fig. 4.2) of 
the beam. The A vector is initialized for all DOF according to 


the following displacements (and slopes), 


v7(&,0) = .1&? 
. (108) 
Vv, (5,0) = 526 
and velocities of displacements and slopes, 
v°(§,0) =0 (109) 
45 (8,0) 
vV,V 
2 
oak 





i. § 


Figure 4.2 Parabolic and Linear Initial 
Conditions Plot Case Study One and Four 
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The prescribed functions of L(t) and its velocity are, 


L(t) =, - .166t 


(110) 


L(t) -.166 


The L(t), function was constructed to permit the beam to be 


drawn half way (10 inches) into the sleeve (hole) in 60 


seconds. Figure (4.3) is a plot of L(t) and velocity, L(t). 


heh 
10. 


60. : 
time 


Figure 4.3 Length and Velocity Function Plot 
Case Study One 
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This problem was solved for a four, eight, and sixteen 
‘element discretization. This was the only investigation of 
grid independence conducted. The results of this 
investigation were discussed in the subsection on spatial grid 


independence. 


C. SECOND CASE STUDY, PARABOLIC FUNCTIONS OF L(t) 

Two studies are conducted where L(t) is prescribed by 
different parabolic functions. If L(t) is defined as a 
parabolic function, its first and second derivatives are non 
zero and significant (See Figures 4.5 and 4.6). Thus, the 
system in Equation (99) is completely defined. In addition to 
the dissimilar functions of L(t), the two cases are distinct 
in their cross-sectional geometries. 

The same initial conditions were imposed for the two 
cases. An initial deformation of the beam in a parabolic 
shape was imposed again as in the first case study and again 
the initial velocities are zero. The A vector of Equation 
(99) was initialized for all DOF according to the following 


displacements and slopes (See Figure 4.4), 


&? 
26 


v"(&, 0) 


(111) 
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and velocities, 


I 
Oo 


v"(§,0) (112) 
¥5(5,0) = 0 , 





T. § 


Figure 4.4 Parabolic and Linear Initial Conditions 
Plots Case Study Two 


1. Parabolic L(t), The Less Stiff Beam 
"Case One" of the second case study is the less 
"stiff" problem. Figure (4.1) shows the cross sectional 
dimensions. These dimensions result in material and geometric 
properties such that p(mass/unit length) = 3.493E-6 1bf: S?/in’, 


and moment of inertia I = 10.17E-6 in’. 








A parabolic function of L(t) is prescribed such that 
the beam is drawn to half its original length in 2.5 seconds, 
reverses direction, and returns to it original length during 
the next 2.5 seconds, for a total of 5 seconds. Figure (4.5) 
is a plot of L(t) and its derivatives, velocity and 


acceleration. The functions are, 


L(t) =L, - 4t + .8t? 


; (113) 
Z(t) = -4+1.6t {0 <t © 5.0} sec. 
L(t) = 1.6 

Ly 

10. | 


. time 








1.6 ; 
0.0 time 
5s 


Figure 4.5 Length, Velocity, and Acceleration Plots 
Case Study Two (Less Stiff Beam) 
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2. Parabolic L(t), The Stiff Beam 

"Case Two" of the second case study is the stiffer of 
the two parabolic L(t) case studies. Figure (4.1) shows the 
cross sectional dimensions of the beam. These dimensions 
result in material and geometric properties such that 
p(mass/unit length) = 6.988E-6 l1bf-S?/in’?, and moment of 
inertia I = 81.38E-6 in‘. 

"Case Two" was started with the same parabolic 
function for L(t) as "Case One". It was here that the 
significance of "stiffness" and computational time came to 
focus. Running the two cases simultaneously as separate jobs 
on different system accounts, clearly demonstrated the 
difference in CPU requirements for the two problems. In fact, 
there was such a disparity in computational effort that it was 
decided to change the course of the stiffer problem ("Case 
Two") such that it’s symmetric, cycle would be complete in 2.1 
seconds vice the 5 seconds of "Case One". The functions of 
L(t) and their derivatives, along with their respective time 


domains are given here, 


L(t) = L, - 4t + .8t? 
L(t) = -4+1.6t {0 $t $1.0} sec. 
L(t) = 1.6 











L(t) = 33.2 - 50.4¢ + 24t? 

L(t) = -50.4 + 48t {1.0 $t $1.1} sec. 
L(t) = 48.0 

L(t) = 5.128 + .64 + .8t? 

L(t) = .64 +1.6t {1.1 S$ t $2.1} sec. 
E(t) = 1.6 


These functions are plotted below in Figure (4.6). 











L 
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. time 
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Figure 4.6 Length, Velocity, and Acceleration Function Plots 
Case Study Two (Stiff Beam) 





3. Discussion of the Parabolic Cases 

The results of these runs provided one of the thought 
provoking questions of the research. The purpose behind the 
parabolic function of L(t) was to observe the beam’s activity 
in the case where it returned to its original length in a 
symmetric, cyclic fashion. The interest was in the question, 
would the expected symmetric behavior of the ‘pull-push’ 
sequence of L(t) be predicted by the code? The results 
clearly showed that symmetry did not occur for the parabolic 
cases (See Figures (4.10) and (4.11)). In fact, the 
deflections had grown considerably during the ‘push’ stage of 
the L(t) cycle. Where did the energy to cause such large 
deflections come from? One possible explanation considered 
was that through the imposition of L(t), energy in the form of 
work had been added to the system. This question was 


addressed in the final case study, wherein the work associated 


with {Faz was tracked. Another possibility, if the work 


cannot account for the increase in displacements, is that the 
results are not correct due to a break down in the numerical 
integration during the latter stage of the ‘push’ stage of the 
cycle. 

Before continuing on to the next set of case studies, 
a discussion of what is a most thought provoking question 
resulting from the research thus far follows. What would 


happen if the beam where drawn totally through the 
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frictionless hole and then pushed back out to its original 
“length? At the end of the ‘pull’ stage of the cycle, the 
entire beam resides motionless in a straight sleeve (the 
frictionless hole) and therefore there is neither 
deformational (strain) energy or vibrating (kinetic) energy. 
Again, energy transfer out of the system as work could account 
for this phenomena. In any case, it may not be possible to 
show this with this numerical model for the following reasons. 
For one thing, as the length of the beam shortens, the shear 
and rotary inertia terms, which were not included in this 
model, become ever increasingly significant and in fact may 
dominate the physical behavior. Secondly, even if the 
physical model could be modified to include these effects, the 
frequencies tend toward infinity as L(t) approaches zero, and 


numerical integration would not be possible. 


D. HARMONIC L(t) PRESCRIBED 

Two studies were conducted simultaneously on two beams 
with the material and geometric properties identical to the 
two beams used in the parabolic L(t) study. In this study, 


L(t) was prescribed as the trigonometric functions of L(t) and 
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its derivatives given here, 


Sage nt 

L(t) 9 cos| "5, 

! - - 7 rt < < 1 
L(t) zig 5n| FS {0 <t $ 3.0} sec. (116) 
si ee as ae t 

ANE ("5 cos| 5] 


Accordingly, L(t) for these cases varied between eight and ten 
inches (Fig. 4.7). The symmetric, cyclic concept was used 


again as it was in the parabolic L(t) prescribed cases. 





tre 


time 


—» time 


Figure 4.7 Length, Velocity, and Acceleration Function Plots 
Case Study Three (Both Beams) 





A sinusoidal function was also prescribed for the initial 
deformation of the 10 inch plexiglass beam. The A vector was 
initialized for all DOF in accordance with the following 


displacements and slopes, 


v'(6,0) = 0.1 - cos| me ) 





2.0 
(117) 
. 3, aed 
vi (6,0) soo SIM FS 
and velocities, 
v°(6,0) =0 (118) 
vi (6,0) = 0 


The initial conditions on deflection and slope are shown 


graphically in Figure (4.8). 





Figure 4.8 Sinusoidal Initial Conditions Plot 
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The previous results of increasing displacements (above 
the initial displacements) observed for the parabolic L(t) 
studies were not obtained in these harmonic L(t) case studies 
(See Figures (4.12) and (4.13)). If the work explanation is 
the correct one in the previous section, then it might be that 
work is associated with parabolic L(t) axial motions, and not 
with harmonic L(t) axial motions. In order to investigate 
this question further, a investigation was undertaken to track 
work for the parabolic L(t) case. This is discussed in the 
next section. The computational effort observations of the 
previous cases where noted again as well. That is, the CPU 
requirement for the "stiffer" problem was greater than for the 
"less stiff" problem, as it had been for the parabolic L(t) 


cases. 


E. CASE STUDY FOUR, TRACKING WORK FOR 2 PARABOLIC L(t) 
A final study was conducted using the "less stiff" beam in 
which a parabolic L(t) was prescribed. The function and its 


derivatives follow and are plotted in Figure (4.9). 


L(t) = 10.0 - 2.666t + .888t? 
L(t) = 2.666 +1.777¢t {0 < t $< 3.0} sec. 
Ett) Si. 777 
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The A vector was initialized for all DOF according to the 
initial deformation of the beam defined by the following 


displacements and slopes, 


v°(&,0) = .1&? 
. (120) 
vi (5,0) = .2€ 
and velocities, 
v"(€,0) = 0 (121) 
v7 (6,0) = 0 


These are the same initial conditions as used in the first 


Case Study and are plotted in Figure (4.2). 





L 

0.0 time 
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9.0 time 
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Figure 4.9 Length, velocity, and Acceleration Plots 
Case Study Four 
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The purpose of this final study was to determine whether 
the increased displacements predicted by the code for the 
parabolic L(t) cases. during the latter stages of the ‘push’ 
stage of the cycle, could be accounted for by work input to 
the system. In addition to tracking work, the moment and 
shear were also tracked. The shear diagrams, shown in Figures 
(4.15) to (4.17), and moment diagrams, shown in Figures (4.18) 
to (4.20), appear to be reasonable. The small values of these 
parameters is due to the values of Young’s modulus, E, and 
moment of inertia, I, used in this study. The product of EI 
for the cases studies here are 1.017 lb. in’, and 8.138 lb. 
Eire ast 

The diagrams for axial force F and work W, shown in Figure 
(4.21) do not appear to be reasonable and therefore are 
suspect. Assuming a one to one relation between L(t) and F(t) 
exists, it is difficult to imagine that such a force would 
produce the smooth parabolic L(t) and vice versa. A tentative 
conclusion therefore is that either F(t) was not coded 
correctly or that there is in fact an instability in the 
numerical integration during the latter stage of the ‘push’ 
cycle of the problem. An effort is presently underway to 
determine if the coding for the calculation of F(t) is 
correct. It should be remarked however that prior to the 
erratic behavior of F(t), which occurs late in the ‘pull-push’ 


cycle, the values of F(t) seemed reasonable. 
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F. FINAL COMMENTS AND RECOMMENDATIONS 

| The results of this initial investigation on the behavior 
of a vibrating beam subject to a prescribed axial motion leads 
to the following conclusions. First and foremost is that the 
implementation of the FEM numerical scheme was accomplished 
with success, although it is not certain that some numerical 
difficulties are not encountered at the later stages of the 
analysis. Further work must be undertaken to resolve whether 
the increase in vibration amplitude which is predicted by the 
code is an actual result of work input to the system or 
whether it is associated with a numerical instability. Prior 
to the investigation of the ‘real’ gun barrel problem, One 
might also investigate whether the omission of axial strain 
energy form the model, which is common whenever bending and 
bar activity coexists, could account for this behavior. If 
so, additional terms for axial strain energy could be included 
in the formulation. 

It is interesting to note that the equation of axial 
motion relates the axial force F(t) not only to the axial 
acceleration L, in accordance with Newton’s law of motion for 
rigid bodies, but also includes additionai terms associated 
with the deformational strain energy of bending, and the 
kineLic energy of beam vibration, at the free end of the beam. 


The former term adds to the acceleration term while the latter 


term decreases it. 
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The practical problem associated with the axial motion of 


a gun barrel due to the recoil action of firing, which 
provided the impetus for this study, was formulated but not 
solved here. An experimental investigation should be 


undertaken to ascertain the accuracy of the numerical model. 
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Figure 4.10 
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Parabolic Axial Motion Transient Response 
Case Study Two (Stiff Beam) 


75 


2.25 


* 
4 xX vr. 





PARABOLIC AXIAL MOTION 


r ’ i : T 1 1 
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 


TIME (SEC) 


Figure 4.11 Parabolic Axial Motion Transient Response 
Case Study Two (Less Stiff Beam) 
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Figure 4.13 Harmonic Axial Motion Transient Response 
Case Study Three (Less Stiff Beam) 
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Figure 4.14 
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Parabolic Axial Motion Transient Response 
Case Study Four 
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Figure 4.15 Shear Plot (0.0 - 1.0 Seconds) 
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SHEAR 


Figure 4.16 
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Shear Plot (1.0 - 2.0 Seconds) 
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Figure 4.18 Moment Plot (0.0 - 1.0 Seconds) 
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Figure 4.19 Moment Plot (1.0 - 2.0 Seconds 
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4.20 Moment Plot (2.0 - 3.0 Seconds) 
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Force and Work Plots 


Figure 4.21 





APPENDIX A 


THE CUBIC SPLINE SHAPE FUNCTIONS 


The beam element is constructed using four shape functions 


(Gir Gezr G3, and q,) which satisfy the following conditions. 
¥ $i q(NP1)=1 , qi(NP1)-0 
\. 
L q,(NP2)=0 , q,(NP2) =0 
’ 2 
‘ ‘hel | 
q,(NP1)=0 , q2 (NP1) =1 
= gq, (NP2) =0 gq, (NP2) =0 
i. 
z 


bs 


J 


q,(NP1)=0 , q3(NP1) =0 


X% gj (NP2)=1 , = gh (NP2) =0 
j 


b+ Q(NP1)=0 , qi (NP1) =0 
oe a q,(NP2)=0 , q(NP2) =1 


where the (’) superscript represents a differentiation with 


Lespect to the spatial variable. 





By sutisfying the four conditions on each q,, the four 


element shape functions can be constructed from a cubic 


equation of the form, 


g,=a,+b,a+c, 0 +d, 


q, is constructed here as an example. 


Gq=a+ba+rce@+d 


1. g,(NP1)=1 > 1 =a, 
2. gi(NP1)=0 > 0=b, 
Y 
gq, = 1+ 6, 0 .+ dd, oF 
\ 


3. g,(NP2)=0 = 0 


1+oc,12+d, 13 


3 


4. gi (NP2)=0 = 0 = 2c, 1,+ 3d, 12 => ¢, = -5da, 1, 
. _ 2 = 3 
substitute (4.) ~ (3.) > d, = = and, c, = -—5 
e le 
thus, 
- 2 2. 
ql az rT a 
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Using the conditions which define shape functions y, q3, 





and, q (listed on the previous page) the other three shape 


functions are obtained, 


gq = 


qd, = 


q, = 





APPENDIX B 


CONSTRUCTION OF THE GLOBAL MATRICES 


The global matrices A, B, C, and, D are constructed from 


the element matrices a*, b*, c*, and, d®* according to 


relationships; 
a-[° o’(97)'a& =Uas 
B= {' & (97) a& = Ub 
(129) 
/ 
c= {* &o(o7)'ab = Ue 
p= {' g(o7) a =Ua 
where, 
/ 
e- i] iT 
a f. q (¢ ) da 
b, =  { q(q7)’ da 
1, 
(130) 


c -Ef, q(q™) ao 


d =f gq (q”) do 


given that 2 is a constant which approximates § transformed 


to the local coordinate @. 
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The transformation of §& is as follows. Referring to 


Figure (B.1), 


s ee ee aes) where, 1, 


i} 
be 


and, £2 = Q? + AN + v2 


Because the transformation of §& to a function of @ results in 


integrals which are difficult to evaluate, it is desireable to 


use an alternative strategy. If we let, 
B= E =n + 2 
MO De 
a numeric value can be assigned to this quantity. Thus, the 


difficult integration is eliminated. Any accuracy lost in the 
approximation will be recovered by additional interations to 


obtain convergence of the FEM solution. 


Figure B.1 
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The final a*, b*, c*, and, d®* matrices follow, 




















12 6 _12 6 
ae Bee ae 
6 4 _6 2 
ap. ts ee 
ae = 
_12 _6 12 _6 
De: ee gat ae 
6 2 _6 4 
ae ee 
1 36- cde G6: ' . 7 
5. TO Bre: 0: 
_1 _22, 1 4, 
bee 20. °° “10° ~36 
6 Ty -2 S6.. | oa 
Si, 40: “Bio 10 
aa bi 1 _2i, 
P20: “30, “Td “as 
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1314, “Hid, 95. 21327 
35 210 70 420 


sis i oe L322 bi 


Zio 105 “420 40 
G1... A3is 13st, id, 
TO 820 TSR TO 
D326 de PhS: 2G 

“420 140 “210 T1065 
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APPENDIX C 


STATIC CANTILEVER BEAM FORTRAN CODE 
AND SAMPLE OUTPUT 


33 3 33 IE 3 IE DEE 3 BE 36 3 3 3 36 3 36 3 38 3 36 3 38 3 IC 3 6 3 IE FE 3 56 JE DEE 3 90 9 C3 EJ 3 3 3 CDE EDC 3 CJ DE JE 9 IE IE 
MARK R. DEVRIES LT USCG 
NAVAL POSTGRADUATE SCHOOL 
SEPTEMBER 1990 
THESIS 
MASTER OF SCIENCE IN MECHANICL ENGINEERING 


TITLE: 
VIBRATION OF A CANTILEVER BEAN 
THAT SLIDES AXIALLY IN A FRICTIONLESS HOLE 


THE FOLLOWING FORTRAN CODE IN A VERIFICATION OF THE FINITE 
ELEMENT FORMULATION FOR THE TRANSIENT PROBLEM TO BE 
PURSUED IN THE NEXT PROGRAMING STEP. 

THE PROGRAM VERIFIES THE [NITE ELEMENT METHOD 

CODE LOGIC ON THREE POSSIBLE STATIC BEAM PROBLEM. 

C1) FIXED END WITH TWO ROLLER SUPPORTS, ONE AT THE CENTER 
AND A SECOND AT THE OPPOSITE END. THIS BEAM IS LOADED BY 
A CONCENTRATED MOMENT AT THE ROLLER SUPPORTED END. 

(2) AND (3) ARE CANTILEVER BEAM PROBLEMS, ONE LOADED BY A 
CONCENTRATED LOAD AT THE FREE END AND THE OTHER LOADED 

BY A CONCENTRATED MOMENT. 

THE PROGRAM IS NOT FLEXIBLE IN THAT IT REQUIRES EDITTING 
AS NOTED IN COMMENT LINES IN THE FOLLOWING SUBROUTINES 
DEPENDING ON WHICH OF THE (3) CASES IS BEING RUN. THE SOLE 
PURPOSE OF THIS PROGRAM IS TO VERIFY THE FEM FORMULATION, 
IT IS NOT INTENDED TO IMPRESS SOLFWARE ENGINEERS. 


(1) SUBROUTINE BC 
(2) SUBROUTINE OUTPUT 
PEP Se LeCE SS SS SSE SOLS SSS SS SS Ste SSS Ee SESS SS ES ESET ESET SSS SS ES St 2 


KOK OK OK OK OF OK OK OK OOK OR OOK OK OK OK OK OK OK OK OOK OK OOK OK OK OK OK OK OK OKOK 


FE HE FE HE OIE HE IE JE FE IE FE SE BE DE IEE DE BE DE BE 36 DE DE JE IE IE 3G HE DE FE IEE 36 3G FE JE IE UE DE BE IE IE DEE DEE IE HE DEE HE De HEE EEK 


% VARIABLE IDENTIFICATION % 
3 3 EE IEE EE 3 IEE 33 EE 3 IE 26 9 EE 3 IE 3 IE HE 3 3 DE KE DEE SE EE 3 BEDE 3 3 EE 3 IE 3 OE DE 3 HE DE IE IE 3 
® NEL - NUMBER OF ELEMENTS 


*% NSNP - RUMBER OF SYSTEM NODAL POINTS 

% NDOF - NUMBER OF DEGREES OF FREEDOM 

* E - MATERIAL MODULUS OF ELASTICITY 

* GI - SECOND MONENT OF THE BEAM CROSS-SECTION AREA 

% BLGTH - BEAM STATIC LENGTH 

® ELE - ELEMENT LENGTH 

* BCM - Se het &PPLIED MOMENT AT FREE OR SIMPLY SUPPORTED 
x EN 

® BCFORC- EXTERNALLY APPLIED FORCE AT FREE OR SIMPLY 

¥ SUPPORTED END 

*% NDETRM - VARIABLE USED IN LOGIC STATEMENT FOR TYPE OF B.C. 
* SLOPE - SLOPE AT FREE END OF CANTILEAVER BEAM 

* DEFLEC - DEFLECTION AT FREE END OF CANTILEAVER BEAM 

*% FACTOR - SCALOR NON-DIMENSIONAL GROUP 

2. 2.2.0,8,.2,2,.5,0,2,2,2,2,5,0,0,2,.2,0,8 5.0,0,6,0,0 5,5,0,0,0,2.5,5,0,0,0.2,2,0,5,2,0,0,5,0,0,0,5,0,0,0.0,2,0.0,0,0,5,5,5 


PARAMETER (N=70) 
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aNANa 


Qa aNN 


OOOO 


20 


25 


26 


27 


28 


29 


40 





DIMENSION ACN,N),FON)D 

DIMENSION WKAREA(6000) 

OPENC10, FILE="/MATRIX OUTPUT’) 

CALL ZERO CA,F) 

CALL DATA CNEL,NSNP,NDOF,€,GI,BLGTH, ELE, BCM, BCFORC) 
ALL MATX (CA,F,ELE,NEL,NDOF) 

CALL BC (A,F,NDOF,NDETRM) 

WRITE €10,%) CCACI,K),K=1,NDOF),J=1,NDOF) 

CALL LEQT2FCA,1,NDOF,N,F,0,HKAKEA, IER) 

CALL OUTPUT (NSNP,ELE,F,BCM,BLGTH,E,GI,NDETRM,NDOF, BCFORC) 
STOP 

END 


ZERO ALL MATRICES 

HE KI IE IE IE DE HE I HEH 5 3 9 IE 3 DE IE 3D EI DE IE DE IE OE 3 DE OE 3 3 IE DE 3 3 9 3 HE EE DE DE EE DE EEE 
SUBROUTINE ZERO (A,F) 

PARAMETER (N=70) 

DIMENSION ACN,ND, FON) 


DO 20 I=1,N 
FCI)=0.0 
DO 10 J=1,N 
ACT, J)=0.0 
CONTINUE 
CONTINUE 
RETURN 


END 


INPUT DATA 
PE ESE SSeS ESSER S ESS ESS SLES ECE SES ESE SETS S SCTE TESS ESTE TSS SS ES 


SUBRCUTINE DATA (NEL,NSNP,NDOF,E,GI,BLGTH,ELE,BCM, BCFORC) 


PRINT *, "ENTER THE NO. OF ELEMENTS TO BE USED IN THE APPROX.' 
READ ®, NEL 

WRITEC6,2€0) NEL 

FORMAT €72X,'hO. OF ELEMENTS IS',15) 

NSNP=NEL+1 

WRITEC6,25) NSNP 

FORMAT(42X,'NO. OF SYSTEM NODAL POINTS IS',15) 
NDOF=2*NSNP 

WRITEC6,26) NDOF 

FORMAT(“2X,'NO. OF D.O.F. IS',15) 

PRINT *®,'THE MODULUS OF ELASTICITY IS?* 

READ ¥*,E 

WRITEC6,27) E 

FORMAT (72X,'MODULUS OF ELASTICITY IS', F10.1) 
PRINT *,'THE SECOND MOMENT OF THE BEAM CROSS-SECTION AREA IS7° 
READ *,GI 

MRITEC6,28) GI 

FORMATC72X,*THE SECOND MOMENT IS',F10.1) 

PRINT *,'THE INITIAL LENGTH OF THE STATIC BEAM IS?" 
READ *, BLGTH 

WRITEC6,29) BLGTH 

FORMATC/2OX, "THE BEAM LENGTH IS',F8.3) 
ELE=1.0/7FLOATCNEL) 

PRINT ®, "ENTER THE VALUE OF THE APPLIED MOMENT! 
READ *,BCM 

WRITE (6,30) BCM 

FORMATC 42x, "MOMENT=',F8.1) 

PRINT *, ‘ENTER THE VALUE OF THE APPLIED FORCE! 
READ *,BCFCRC 

WRITE (6,40) BCFORC 

FORMAT(42X,*FORCE=',F8.1} 


RETURN 
END 


FILL LARGE A 
PPS ESL See CT eC ese sess ese ces ese ress srs ese teres esses eres ess: 


SUBROUTINE MATX (A.F,ELE,NEL,RDOF) 
PARAMETER CH=70) 
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‘ DIMENSION ACN,N),FCND 
. CALCULATE LITTLE A MATRIX 


A11=12.0/CELE¥X3.0) 
Al2=6.0/(ELEX*2.0) 
A13=(-1.0)Al11 
A14=A12 

A21=A12 
A22=4.0/7ELE 
A23=(-1.0)%A12 
424=2.0/7ELE 
A31=A13 

A352=A23 

A33=Al] 

A34=A23 e 
AG1=A14 

A42=A24 

AG3=A34 

A44=A22 


FILL LARGE A MATRIX 
=2*NEL-1 
DO 10 I=l,L,2 


ACI,I)=A 
ACI,1I+1) 
ACI,I+2) 
ACI,I+33 


ACT4#1,1)= 
ACI+1],1+1 
ACI+1,1+2 
ACI+1,1+3 


a a 9000 


H+ 


+] )+A12 


( 


+1,1)+A21 
C141, 14194422 
24 


wevyp rPpPpr 


C 
10 CONTINUE 


RETURN 
END 


HEHE IE HK IE EE 33 3 HE BE 3 HC BC DE 3 IE 3 BE IE 3 BEE EE 3 IE BE BE IE 33 IE BE 3K IE 9 DE IE 3 CE DF 3 3 ME BE I HE 3 3 HE 
THE FOLLOWING SUBROUTINE ALTERS THE GLOBAL MATRICES TO 
IMPOSE THE BOUNDARY CONDITIONS 


SUBROUTINE BC (A,F,NDOF,NDETRM) 

PARAMETER(N=70) 

DIMENSION ACN,N),FON) 

PRINT *, ENTER 1 FOR THE OVER DETERMINANT CASE OR! 
PRINT *%,'2 FOR THE FREE END CASE.‘ 

READ *, NDETRM 

IF (NDETRM .NE. 1) GOTO 20 


AMEND A TG ACCOUNT FOR BOUNDARY CONDITIONS 


CHANGE FIRST AND SECOND ROHS TO ACCOUNT FOR THE ESSENTIAL 
BOUNDARY CONDITIONS AT THE FIXED END. 


ANANAAN 


ANANQANNAAND 


THE J-TH EQUATION IS THE EQUATION DESCRIBING DEFLECTION AT 
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ANAANNNAN 


oO 


an QANQAANNAG NAOANAAN AN|ANND 


NANNAAMANANANA 


10 


20 


30 


40 


THE LOCATION OF THE CENTER SUPPORT. IT IS REPLACED BY THE 
ESSENTIAL BOUNDARY CONDITION ON DEFLECTION. 


THE CNDOF-1)TH EQUATION IS THE EQUATION DESCRIBING THE 
DEFLECTION AT THE ROLLER SUPPORTED END. THIS EQUATION 
IS REPLACED 2Y THE ESSENTIAL B.C. ON DEFLECTION. 


J=NDOF/Z2 
DO 10 I= 


»HDOF-1)=1.0 


¥RXXPOINT LOADKRR 
THIS LINE IS ACTIVATED FOR THE CANTILEVER BEAM LOADED 
BY A CONCENTRATED LOAD CASE 


FCNDOF-1)=1.0 


¥*XXXPOINT MOMENT ¥*x 
THIS LINE MUST BE ACTIVATED FOR BOTH THE OVER DETERMINATE CASE 
AND THE CANTILEVER BEAM WITH A CONCENTRATED MOMENT CASE 


FOCNDOF)=1.0 


RETURN 
END 


FORMULATE OUTPUT 

PERE E Cee eee Sete S i SESE SSS SESS ES SSE ST ES ES ES SETS SESE TESS SSS 
THIS SUBROUTINE CALCULATES THE EXACT SOLUTION FOR THE 
Coe BEAM CASES AS WELL AS PRINTS THE OUTPUT OF ALL 
CASES. 


SUBROUTINE OUTPUT CNSNP,ELE,F,BCM, BLGTH,E,GI,NDETRM,NDOF, BCFORC) 
PARAMETERCN=70) 
DIMENSION FCN) 


X¥¥POINT FORCE RX 
FACTGR=BCFORCX( BLGTH¥*®2.0)/CExXGI) 
DO 5 I=1,NDOF-1,2 
FCI) =FCI)*¥BLGTH¥FACTOR 
FCI+1)=FCI+1) FACTOR 
CONTINUE 


¥¥¥POINT MOMENT *%% 
FACTOR=BLGTH¥BCM/ (CEXGI ) 

DO 5 I=1,NDOF-1,2 
FCI)=FCI)*XBLGTH*¥FACTOR 
FCI+1)=FCI+1) FACTOR 

CONTINUE 


WRITEC, 30) 
WRITECX, 40) 
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iE ESSSSSSSSSSSSSS SIS SS aC 


agNgnNg 


ANQNAMANN 


10 


15 


15 


IF (NDETRM .NE. 1) GOTO 15 
WRITEC%, 50) 
Jel 
XLOC=0.0 
DO 10 I=1,NSNP,1 
HWRITEC*,20) XLOC, FCJ), FCS+1) 
XLOC=XLOC+ELE 
JeJ+2 
CONTINUE 
GO TO 80 


HNXNXPOINT LOAD HX 
ACTIVATE FOR THE CANTILEVER-CONCENTRATE FORCE CASE 


SLOPE=(BCFORCK( BLGTHX¥®2.0))/(2.OXEXGI) 7 
DEFLEC=( BCFORCX( BLGTH*X3.0))/(3.0XEXGI) 


*X*¥XPOINT MOMENT XX XX 
ACTIVATE FOR THE CANTILEVER BEAM CONCENTRATED MOMENT CASE 


SLOPE=(BCMXBLGTH)/CEXGI) 
DEFLEC=( BCMX( BLGTHXX2.0))702.0XEXGI) 


WRITEC#,60) 
WRITEC*,65) 
WRITEC*,70) DEFLEC,FCNDOF-1),SLOPE, FCNDOF) 


FORMAT (2X,F8.3,3X,E£12.4,3X,E12.49) 

FORMATC2X,° X-LOCATION', 2X, "DEFLECTION', 3X, *SLOPE' ) 
FORMATCIX,! DEFLECTION AT B ees 

FORMATC6X, "EXACT, 11X, "FEM", 11X, "EXACT", 11X,' FEM") 
FORMATCIX,E12.4,3X,E12.9,3X,E12.4,3xX,E12.4) 

RETURN 

END 


SLOPE AT B') 





NO. OF ELEMENTS IS 8 
NO. OF SYSTEM NODAL POINTS IS 9 

: NO. OF D.O.F. IS 18 
MODULUS OF ELASTICITY IS 30000000.0 
THE SECOND MOMENT IS 100.0 
THE BEAM LENGTH IS 100.000 

. ENTER THE VALUE OF THE APPLIED MOMENT 
MOMENT = 0.0 
ENTER THE VALUE OF THE APPLIED FORCE 
FORCE= 1000.0 
ENTER 1 FOR THE OVER DETERMINANT CASE OR 
2 FOR THE FREE END CASE. 
? 


2 


DEFLECTION AT B SLOPE AT B 
ACT FEM 


EX EXACT 
0.1111E+00 O.1111E+00 0.1667E-02 
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FEM 
0.1667E-02 


APPENDIX D 


INTEGRATION BY PARTS 


The following is the detail of the integration by parts on 


the first term of Equation (54), 


1 
J, 2 (07 hee B 8 (137) 


The first integration results in, 


1 
f, 2 (9 ee d§ d= 
(138) 


1 
o(97).. 3 - fo, (07), d& 8 


A second integration performed on the integral in Equation (2) 


gives, 
1 
- {° 0 (0 he a6 8 = 


; (139) 
- 9, (07), 813+ f° @, (07), ab 8 
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Combining Equations (1), (2), and (3) gives the symmetric 


operator and boundary terms below, 


. I 2 (97 ee dé 6 = 


i. ff (140) 
. lo (9° has 5 - o (9 he : | : i ats (o he ae 
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APPENDIX E 


TRANSIENT BEHAVIOR OF A CANTILEVER BEAM FORTRAN CODE 7 


HEE IE 3 3 HE IE 3 HE <3 3 HE IE 3 HE DE 3 HE DE BE IE DE 3 HE BE BE 3 DE BE 3 HE BE 3 HE 3 3 IE HE DE IC HE BE DE IE BE IE IE 3 HE DE DE CHE DE IE WF HEHE 
MARK R DEVRIES LT USCG 
NAVAL POSTGRAQUATE SCHOOL 
SEPTEMBER 1990 


THESIS 
MASTER OF SCIENCE IN MECHANICAL ENGINEERING 


THAT SLIDES AXIALLY IN A FRICTIONLESS HOLE 


x 

% 

x 

5 3 

x 

x 

% TITLE: 

* VIBRATION OF A CANTILEVER BEAM 

x 

% 

% THE FOLLOWING FORTRAN CODE UTILIZES THE FINITE ELEMENT 

% METHOD AND AN IMSL PACKAGE INTEGRATION SUBROUTINE DIVPAG 

* TO SOLVE THE ABOVE PROBLEM. THE PROGRAM IS HIRITTEN WITH 

%® NUMCROUS COMMENT LINES WHICH EXPLAIN THE CODING. 

BE HE IE IE IE I HK HE HE IE HE HE HE HE HE HK HE HE HE HE HE HE HE HE HE IE IE HE HE HE HE IE HE DE HE HEHE HE HE HE HEE HE HEH HEE HEE HE HE EEE 
eee SS SSS SSS SESS SSS SSS SSS SESS SSS SS SS ESS SOS S SSS SESE SSS SSE SS SS | 
¥ VARIABLE IDENTIFICATION * 
3 HE IE IE HE IK FE IE BE FE IE 3K FE BE FE SE DE HE 3K FE IEE OE HE SE HE HE FE IE HE IE HE FE HE DE HE IE HE IE HE HEE HE HE KKK HK KK EK HHH KKK 
* NEL - NUMBER OF ELEMENTS 

*% NSNP - NUMBER OF SYSTEM NODAL POINTS 

% NDOF - NUMBER OF DEGREES OF FREEDOM 

% N,NN - DIMENSIONS OF MATRICES AS SPECIFIED IN DIMENSION 

¥ STATEMENTS : 
% €& - MATERIAL MODULUS OF ELASTICITY 

* GI - SECOND MOMENT OF THE BEAM CROSS-SECTION AREA 

* ELE - ELEMENT LENGTH 

% ALPHA - LOCATION OF ELEMENT LEFT GNP » 
* PSIAVE -_ ESTIMATE OF PSI 

* PSISQ - ESTIMATE OF PSI SQUARED 

* TEND ~ VALUE GF TIME AT WHICH THE SOLUTION IS DESIRED 

®% NEQ - NO. OF FIRST ORDER DIFFERENTIAL EQUATIONS 

*% TIME - INDEPENDENT TIME VARIABLE 

* DELTME - TOTAL TIME INCREMENT FOR ONE INTEGRATION STEP 

* BETA - CONSTANT DETERMINED BY BEAM MATERIAL PROPERTIES ONLY 
®% RATE - LENGTH CHANGE PER UNIT TIME 

®* EXEE - THE GLOBAL NONDIMENSIONAL AXIS, THAT IS, (X/L) 

* DELTA - THE VECTOR OF NONDIMENSIONAL NODE DEFLECTIONS 

x AND SLOPES. MUST BE MULTIPLIED BY L(T) FOR 

* ACTUAL DEFLECTIONS. SLOPES REMAIN THE SAME. 

x 
x 


IE IE 3 3 3 HE 9 DE 3 3 3 3 9K 3 3 3 IE 3 3 3 3 9 3 3 3 3 3 3 3 3 CE 3 0 3 0 3 3 0 9 9 EE DE 3 IEE 
INCLUDE "COMMON FORTRAN’ 
DIMENSION CELTACNN), PARAMCNPARAM) ,WKSCNN) 
DIMENSION YPRIMECNN) 
COMMON /“WORKSP/ RWKSP 


REAL RWKSP(6608) 


EXTERNAL FCN 
EXTERNAL FCNJ 








(olere) 





¥X 
RR 
*X 
KR 
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OPENC9, FILE="/MARKI INPUT") 
OPENC10, FILE="/MARK3B OUTPUT’) 
OPENC11, FILE="/DATA] INPUT) 
OPEN(C12, FILE='/DATA2 INPUT") 


CALL IWKIN(G6608) 

CALL ZERO (DELTA, PARAM) 
CALL DATA 

CALL MATX 


PI = 3.141592654 


DEFINITION OF PARAMETERS REQUIRED BY IMSL MATH LIBRARY ROUTINE 
DIVPAG 


IDO=1 

NEQ=2*NDOF 
TIME=TSTART 
TOL=1.0E-4 
PARAM(4)=2000000 


PARAMC12) IS 1 FOR ADAMS METHOD AND 2 FOR GEAR (STIFF) METHOD 


PARAM(12)=2 
PARAM(13)=2 
PARANC19)=1 
PARAN(20)=NN 


INITIALIZE THE DEPENDENT VARIABLE ARRAY DELTACNEQ). 


CAUTION: THE NONDIMENSIONAL VSTAR IS CONSTRUCTED HERE. 

TO OBTAIN THE ACTUAL INITIAL DISPLACEMENT CONFIGURATION, V, 
SUBSTITUTE THE NONDIMENSIONAL COORDINATE AXIS EXEE IN THE 
EXPRESSIONS BELOW BY (X/ZLINT) AND REPLACE VSTAR*ZLINT BY V 


IF CISTART.EQ.0) THEN 
IC=NDOF-1 
EXEE = 0.0 
MRITEC1O,%) "THE INITIAL TIME PRIOR TO INTEGRATION = ',TIME 
WRITEC1]O,%) "THE INITIAL DELTA VECTOR IS 


DO 10 I=1,I1C,2 

PIOV2 = PI/2. 

DELTACI) = 0.1 - 0.1*COSCPIOV2XEXEE) 
DELTACI+1) = O.1¥PIOV2XSINCPIOV2XEXEE) 
DELTACI)=0. 1 XC EXEEX%2.0) 

DELTACI+1) = 0.2xEXEE 


WRITE €10,%) DELTACT) 
VIRITE €10,%) DELTACI+1) 
EXEE = EXEE + ELE 


u 


CONTINUE 
ELSE 
READ(11,¥) TSTART,ZL 
WRITECLO,%) "RESTART TIME = ',TSTART, "WITH LENGTH ',2L 
READ(12,%) CDELTACII), II=1,NEQ) 
WRITECLO,%) ‘INITIAL DELTA IN RESTART FOLLONS!' 
WRITEC1O,*) CDELTACJJ),JJ=1,NEQ) 
TIME=TSTART 
END IF 
IDO=1] 


CALL FORMCHEQ, TIME) 


pO 10900 IEQ=1,N 
YPRIMECTEQ)= 
DQ 900 IC=l, 


EQ 
0. 
NE 


0 
Q 








Cc 
c 


900 
1000 


30 





YPRIMECIEQ) =YPRIMECIEQ)+HCIEQ, IC)XDELTACIC) 
CONTINUE 
CONTINUE 


WRITEC10,%) ‘INITIAL YPRIME VECTOR PRIOR TO ENTRY TO DIVPAG' 
WRITE (10, oe a , TQ=1,NEQ) 


WRITECIO, 
WRITEC1O, baa 
WRITECLO, ‘Il) INTEGRATION LOOP, TIME, LENGTH AND DELTA FOLLOWS’ 


WORK = 0.0 

FL = RHOXZLINT®2.*ACC 

VEE = 6.*(DELTACNDOF-3)-DELTACNDOF-1))/ELE2 
& + 2.¥*CDELTACNDOF-2)+2.¥*DELTACNDOFIIZELE 

F2 = CEXGID¥CVEEXX2)/7CZLINTX¥2) 

F3 = RHOXC-RATEXDELTACNDOF) + ZLINTXDELTACNEQ-1 ) ) x¥2 
FHEW = Fl + 0.5*C(F2 - F3) 

ZL = ZLINT 


DO 30 JEND=1,NSTEP 

FOLD = FNEW 

ZLOLD = 2L 
TEND=TSTART+DELTMEXFLOATCIEND)/FLOATCNSTEP) 
IFCTEND.GT.3.0) GO TO 35 

CALL DTIMECIHOUR, MINUTE, ISEC) 

IFCIHOUR .LT. 18 .AND. IHOUR .GE. 7) GO TO 35 
CALL DIVPAG CIDO,NEQ,FCN,FCNJ,G, TIME, TEND, TOL, PARAM, DELTA) 
IF (MODCTEND,1).EQ.0) THEN 

ZL = ZLINT - RATEXTIME + ACCKCTIMEX#2) 

ZL2 = ZL¥¥2 

ZLDOT = -RATE + 2.¥*ACCXTIME 

ZL DDOT = 2.x*ACC 

ZL = 9. + 1.% COS(€ PIXTIME/1.5) 


WRITEC1O,*) ' ' 

WRITECLO,*) ¢ ! 

WRITECIO,*) "TIME = ', TIME, ‘LENGTH = ', 2L 

HRITEC1O,*) *% 3 

HWRITEC1O,%3 "DELTA FOLLOWS! 

WRITECLO,%) CDELTACIQ), 1Q=1,NEQ) c 
Z2MOM = CEXGI/ZL) RC C6. Fate te DELTAC1)) 

& -(2.7ELE)*CDELTAC4)+2.xDELTAC2) 3 


SHEAR = CEXGI/ZL2)¥¢ (12,7 ELE3) €¢ DELTACL)=DELTAC3)) : 
& +(6./ELE2) XCDELTAC2)+DELTAC4G) )) 

WRITEC1O,*) ' ' 

ere "MOMENT = ', ZMOM, "SHEAR = ', SHEAR 

END 


Fl = RHOXZLINTX2. ACC 

VEE = 6.*XC(DELTACNDOF-3)-DELTACNDOF-1))/ELE2 
& + 2.*CDELTACNDOF-2)+2.*DELTACNDOF))/ELE 
= CEXGI)*XCVEEXX2)7ZL2 
F3 = RHOXC-ZLDOT¥DELTACNDOF) + ZL¥DELTACNEQ-1) ) x2 
FNEW = Fl + G0.5x*CF2 - F3) 


DELWORK = O.5x®CFNEW + FOLD) *(Z2L - ZLOLD?) 


WORK = WORK + DELWORK 
WRITEC1O,*) 'ZLDOT ', ZLDOT, *ZLDDOT = ',ZLDDOT 


HIRITEC1O,*) ‘OLD F = ", FOLD, "NEW F = ', FNEW 

WRITEC10,%) "WORK = *, WORK 

CALL FORM(NEQ, TIME) 

CCNTINUE 

CONTINUE \ 
Ip0=3 

CALL DIVPAG (IDO,NEQ,FCN,FCNJ,G, TIME, TEND, TOL,PARAM, DELTA) 


STOP 
END 
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35 
40 


ZERO ALL MATRICES 
HE HIE IE 5K HE DEH HEHE DE HE IK He IE DEE HE BE HE EH HE BE 3 HE DE 3 BE BE 3 HE I DE IE 3 3 HE BE DE HE DE HE DO HC 3 BE HE HE IE HEHE 


SUBROUTINE ZERO (DELTA,PARAM) 


INCLUDE 'COMMON FORTRAN 
DIMENSION DELTACNN), PARAMCNPARAM) 


moecuunwuu 
Nee ee eee es tt 
wun wns 
Q2oo000Z 
229000 


DO 40 I=1,NN 
DELTACI)= 
BO 


CON 
CONTINUE 
DO 50 I=1,NPARAM 

PARAM(I)=0.0 
CONTINUE 


RETURN 
END 


INPUT DATA 

HIE I HE IEE I HE IEE HEH IE HE DEH HE HE I HE HE IE IE 3 3 IE IE IE DE 3 HE JE IE DE 3 HE 3 9 39 HEE DI HIE 
SUBROUTINE DATA 

INCLUDE "COMMON FORTRAN 

READ (9,%) NEL,E,GI,RHO,NSTEP, ISTART, TSTART,ZLINT, DELTME,RATE,ACC 
WRITEC6,%) "THE NUMBER OF ELEMENTS IS ', NEL 


NSNP=HEL+1 
NDOF=2*NSNP 


WRITEC6, ®) 


HRITEC6, ¥) 
WRITEC6, *) 


"THE 


"THE 
‘THE 
"THE 
*THE 


NUMBER OF SYSTEM NODAL POINTS IS ', NSNP 
NUMBER OF DEGREES OF FREEDOM IS ', NDOF 
MODULUS OF ELASTICITY IS ', & 

WRITEC6, *®) MOMENT OF INERTIA IS ', GI 
WRITECS, ®) MASS PER UNIT LENGTH IS '*, 


ELE=1.0/FLOATCNEL) 
ELE2 = ELEx*x2? 

ELE3 = ELEx*3 
BETA=RHO/CEXGI) 


RHO 


WRITEC6,*) "THE VALUE OF BETA IS * ,BETA 

WRITEC6,%) "THE NUMBER OF INTEGRATION STEPS IS ', NSTEP 
WRITEC6,%) "ISTART I3 1 FSR RESTART; HERE IT IS ', ISTART 
WRITEC6,*) "THE INITIAL LENGTH IS ', ZLINT 

WRITEC6,%) "RATE OF AXIAL MOTION IS * , RATE 

RETURN 

END 


FILL LARGE A 
PE TESST CSCL SESS SS SES ES ESET SSCS TCC Se esse ee Tee Sse ese eee e sees | 


SUBROUTINE MATX 
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Cc 


c 


c 
Cc 


agg 


INCLUDE ‘COMMON FORTRAN' 
CALCULATE LITTLE A MATRIX 


A11=12.0/CELEX*3.0) 
Al12=6 .0/CELE**2.0) 
Al13=(-1.0)8A11 
ALG=AL2 

A21=A12 
A22=4.0/7ELE 
AZ235=(€-1.0)*A12 
A2G=2.0/7ELE 
AZ1=A1L3 

A32=A23 

A33=A11 

A34=A23 

AG1=A1S 

AG42=A24 

AG3=A34 

AGG=A22 


CALCULATE THE ELEMENTAL B MATRIX 


Bll=(€-6.0)/7(5.0XELE) 
Bl2=-1.1 
B13=¢€-1.0)%B1l 
B1l4=-.1 


B22=(-2.Q0ELE)/415.0 
B23=.1 

B24=ELE730.0 
B31=B13 


B44=B22 


CALCULATE THE ELEMENTAL C MATRIX 


Cll=-.5 
C1l2=ELE/10.0 
C13=.5 
C1G=(-.1)*ELE 
C21=C14 
C22=0.0 
C23=C12 
C24G=CELEXX2)7(-60.0) 
C31=Cll 
C32=C21 
C33=C13 
C34=C12 
C41=Cl2 
C42=20-1,.0)¥*C24 
C43=C32 
644=C22 


CALCULATE THE ELEMENTAL D MATRIX 


DI1=13.0*ELE755.0 
D12=(€11.07210.0)xCELEX*2.0) 
D1S=9.Q0¥ELE/770.0 
BPIG=(-13.0/420.0)*( ELEX*2) 
D2e1=d12 

225 CELE*R*S 0970105.) 
D25=(13.O0*ELEx*2)7420. 
C29=~CELEXX3.0)7(140.0) 


7v- 


vain 


13 
32=D23 


o 
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D33=Dll 
D34=(-11 .OXELEX*2)/210. 
D41=D14 
D42=D24 
D43=D34 
D44=D22 


c FILL THE GLOBAL A,B,C, AND D MATRICES 
LNEL = 2¥*NEL-1 
DO 10 I=l, LNEL, 2 
ALPHA=0.0 


PSIAVE=ALPHA+CELE/2.0) 
PSISQ=PSIAVEX*2 


ACI, I)=ACI,1I)+Al1 
ACI, I+t1l)=ACI,I+1)+Al12 
ACT, I+2)=A13 

ACI, I+3)=A14 


ACT4+1,1)=ACI+1,1)+A21 
ACI4+1,I+1)=ACI+1,1+1)+A22 
ACI+1,1+2)=A23 
ACI+1,1+3)=A24 


ACI+2,1)=A31 

ACI+2,1+1)=A32 
ACT#+2,1+2)=A33 
ACI+2,1+3)=A34 


ACT4+3,1)=AG41 

ACI+3,1+1)=AG2 
ACT+3,14+2)=A43 

ACT+3,1+3)=A44 

B11 ¥PSISQ) 
oe eae 


LNA NEN 
++ +r 


Hh Hee WNre u 


yew ee 


+1,0)+C(B21¥PSISQ) 
CI+1,1+1)+¢CBe2*PSisQ) 
B23*PSISQ 

=B2GXPSISQ 


YEBSL¥PSISQ 


»1+2)=B33*PSISQ 
2,1+3)=B34*PSISQ 


+3, 1)=BG1¥PSISQ 
+3,1+1L)=B42*PSISQ 
+3,1+2)=B43*PSISQ i 
+3, 1+3)=BGG*PSI5Q 


I)=CCI,1I+0C1L1¥PSTAVE) 

I+] =CCI,1T+1)2+CC1l2O*PSIAVE) 
I+2)=CL3¥*PSIAVE 
I+ 
1 


I 
I 
3 
4 
I 
B 


++ tet 
NNR Ree ete 
~~ yee 


nannn 


~ 


+ 
rm 
~ 
Les] 
wG 
N 
* 
vu 
“n 
H 
172) 
w%) 


+ 
+ 
+ 
+ 


RANA 


AON 
_ ot et et ot toon Moen Aeon oe I eter tt et he | oe oe oe oe) 


+ 


, 
, 
, 
? 


3)=C1G¥PSIAVE 


io) ANNA wow mnwow wouwww wvwww 
LX ONIONS 


¢ »T)=CCI+1, 1)+0C21¥PSIAVE) 
CCI+l, T+L =CCI+1, 1+1)+0C22xPSrAve) 
CCI+1,[+2)= C23xPSIAVE 

C141, 1+3)=C24XPSIAVE 


CC1+2,1)=C31¥PSTAVE 


CCL+2, 1+] )=C32*PSIAVE 
CCI+2,1+2)=C33*PSIAVE 


107 





CCI+2,1+3)=C34*PSIAVE 


CCI+3, 1 =CGIXPSIAVE 

CCI+3,1+1)=C42"PSIAVE 
CCI+3,1+2)=C43¥PSIAVE 
CCI+3,1+3)=C4GRPSIAVE 


DCI,I)=DCI,1)3+D11 
DCI, 1+1)=DCI,1+1)+D12 
DCI,1+2)=D13 
D(I,1+3)=D14 


DCI+1,1I)=DCI+1,1)+021 
DCI+l, I+1)=DCI+1,1+1)+BD22 
Di I+1,1+2)=D23 
D(I+1,1+33=D24 


D(I+2,1)=D31 

DC I+2,1+1)=D32 
DC I+2,1+2)3=D53 
DC1I+2,1+3)=D34 


DCI+3,1)=D41 

D(1+3,1+1)=D42 
DC1+3,1+2)=D43 
DC 1I+3,1+3)=D44 


ALPHA=ALPHA+tELE 
1Q) CONTINUE 


RETURN 
END 


9333333336 3333093039933 OOOERE 
SUBROUTINE FORMCNEQ, TIME) 
INCLUDE ‘COMMON FORTRAN® 


ZU = ZLINT ~ RATEXTIME + ACCHCTIMERX2) 
ZLDOT = -RATE + 2.*ACCHTIME 

ZLDDOT = 2.¥*ACC 

FI = 3.141592654 

PARAM = PI/1.5 

ZL= 9. + 1.% COSC PARAM*TIME) 

ZLDOT = -CPARAM)*SINCPARAM*®TIME) 
ZLCDOT = -CPARAM*®¥2)*COSCPARAMXT IME) 
ACGEFF=(-1.0)7¢00CZL**G)*BETA) 

ECOEFFS -CCZLDOT/ZL)**2) 

CCGEFF= -.2.xCZLDOT/“ZLIXX2) + ZLDDOT 


PO 20 I=1,NDOF 
CO i5 J=1,NDOF 
RCI.J)= ACT, J)*ACOEFF + BCI,J)*BCOEFF + CCL,J)*CCOEFF 
CCCI,J)= -2.*BETAR(ZLDOTZZL)xCCI, J) 
CONTINUE 
CONTINUE 


aan 


ow 


Ne 


Cc REDUCE SYSTEM OF EQUATIGNS FROM SECOND TO FIRST O.D.E. 
DO 100 I=1,NDOF 
6(I,1)=1.0 
100 CONTINUE 


DO 300 I=1,NDOF 


K=I+NDOF 

DG 200 J=1,NDOF 
LL=J+NDOF 
GCK,LLI=DCI,J) 


co CONTINUE 
00 CONTINUE 
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DO 400 I=1,NDOF 
IT=1+NDOF 
HCI, II)=1.0 
400 CONTINUE 


BDO 600 [=i,NDOF 
K=I+NDOF 
DO 500 J=1,NDOF 
HOCK, JJ=RCI,J) 
CONTINUE 
CONTINUE 


DO 700 I=1,NDOF 
ITI=I+NDOF 
DO 650 J=1,D0F 
JJ JItNDOF 
HCII,JJ) = CCCI,J) 
650 CONTINUE 
700 CONTINUE 


IMPOSE FIXED END BOUNDARY CONDITIONS 


Do 800 J=1,2 
DO 750 I=1,NEQ 
GCJ+HDOF,1)=0.0 
HCJ,1)2=0.0 
HCJ+NDOF,1)=0.0 
750 CONTINUE 
800 CONTINUE 


nut 
Qo 
aa 


GCNDOF+1,NDOF+1)=1. 
GCNDOF+2,NDOF+2)=1. 


RETUPN 
END 


0 
0 


3 36 33 HE HE IE IE I 3 3 IE IE DE BE 3 3 3 3 3 3 3 EE 3 3 3 EE 3 3 3 3 3 9 3 3 3 EE 3 3 EC 3 3 3 3 3 3 3 3 EE 3 3 HEE 
SUBROUTINE FCN (NEQ,TIME, DELTA, YPRIME) 


INCLUDE "COMMON FORTRAN? 
DIMENSION YPRIMECNEQ), DELTACNEQ) 


REAL L 

DELTAC1) = 0.0 
DELTAC2) = 0.0 
DELTACHDOF+1) = 0.0 
DELTACNDOF+2) = 0.0 


FORM YPRIME 


DO 1000 ITEQ=1,NEQ 
YPRIMECTEQ)=0.0 
DO 900 IC=1,NEQ 
YPRIMECIEQ)=YPRIMECIEQ)+HCIEQ, IC)XDELTACIC) 
$00 CONTINUE 
1000 CONTINUE 


RETURN 
END 


FUNCTION FCNJCNEQ, TIME. DELTA, PD) 
REAL “IME, DELTACNEQ), PD(%) 
FCHJ=0.0 

RETURN 

END 
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